Consequences of imperfect mixing in the Gray-Scott model 
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We study an autocatalytic reaction-difTusion scheme, the Gray-Scott model, when the mixing 
processes do not homogenize the reactants. Starting from the master equation, we derive the 
resulting coupled, nonlinear, stochastic partial diflterential equations that rigorously include the 
spatio-temporal fluctuations resulting from the interplay between the reaction and mixing processes. 
The fields are complex and depend on correlated complex noise terms. We implement a novel method 
to solve for these complex fields numerically and extract accurate information about the system 
evolution and stationary states under different mixing regimes. Through this example, we show how 
the reaction induced fluctuations interact with the temporal nonlinearities leading to results that 
differ significantly from the mean-field (perfectly mixed) approach. This procedure can be applied 
to an arbitrary non-linear reaction diffusion scheme. 

PACS numbers: 05.10.Cc, 82.40.Bj, 05.40.Ca, 89.75.Kd 
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Reaction diffusion networks are relevant to a broad 
scope of chemical, biological and nonlinear optical sys- 
tems where the agents must diffuse about before they 
meet and react To date, much of the theoretical 
work devoted to the analysis of reaction network dynam- 
ics has been based on the assumption that the system is 
"well mixed" , i. e. that the concentration of each species 
is always uniform throughout the system. This regime is 
then solved using mean-field approximations which may 
be adequate under perfect stirring conditions or cquiv- 
alently when the diffusion rates are very high. But if 
the reaction rates are nonlinear, which is the case of the 
ubiquitous autocatalytic and catalytic networks, any de- 
parture from the assumption of perfect mixing can have 
major dynamical consequences. There are numerous ex- 
perimental evidences on the outcome of catalytic pro- 
cesses under imperfect mixing Q. For instance, it has 
been experimentally shown that imperfect mixing in au- 
tocatalytic systems can lead to chemical reactions that 
occur at seemingly random intervals crystallization 
processes that yield sometimes all left-handed and some- 
times all right-handed crystals Q and reactions in which 
changing the rate at which a solution is stirred can cause 
a transition from a stationary, time-independent state to 
one of periodic or even chaotic oscillation d. These phe- 
nomena have been studied primarily in inorganic chemi- 
cal media, but there are implications also for relevant bi- 
ological systems: hypercycle networks, viral quasispecies 
dynamics, and ecosystems, etc. [1, 0, There is a 
wealth of theoretical works related to the effects of fluc- 
tuations on these systems that intuitively show the rel- 
evance that density fluctuations may have on stability, 
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pattern formation and the effects of dimensionality 
However, in these approaches, the fluctuations are sim- 
ply added ad-hoc to the mean-field description and the 
noise strength and correlations are chosen arbitrarily. 

In this Brief Report we address this issue and work 
through a general-purpose, mathematical and numeri- 
cal procedure to study the mixing process in reaction- 
diffusion problems. The mixing process in nonlinear re- 
action systems has been described mostly for simple re- 
action schemes [ic|] but has never been fully exploited 
to solve for the Langevin-type equation, since the fluc- 
tuations in the general case, turn out to be complex. 
Our approach is based on (1) the application of a stan- 
dard, fleld-theoretic method to characterize the reaction 
induced fluctuations and derive a Langevin-type descrip- 
tion of the underlying molecular dynamics, (2) the non- 
dimensionalization of the problem to characterize the de- 
pendence on the agents' diffusivities, reaction rates and 
the spatial dimension and (3) the separation of the real 
and imaginary parts of the noise to solve numerically for 
the real and imaginary parts of the stochastic flelds. As 
an example of this consistent description, we treat the 
Gray-Scott (GS) model in two spatial dimensions [ll| 
and vary the diffusion rates to rigorously explore differ- 
ent mixing regimes. We choose this model to facilitate 
visualizing the differences between perfect and imperfect 
mixing as the GS model is known to exhibit spatial pat- 
tern formation and to be very sensitive to multiplicative 
noise [l^- This procedure can be applied to an arbi- 
trary reaction scheme for any spatial dimension. The 
method presented here allows one to solve for the full set 
of complex Langevin equations and to extract information 
about the spatio-temporal time evolution of the system 
and its stationary states under different mixing regimes. 
We show numerically that only the noise averaged real 

The reactants. We compare our results with the mean- 
field case, which is recovered only in the infinite diffusion 
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limit. 

The GS model, corresponds to the following chemical 
reactions: 

U + 2V ^ W, 

Mo 



u. 



(1) 



The concentrations of the chemical species V and U are 
functions of d-dimensional space x and time t. A is the 
reaction rate, p and q are inert products, is the decay 
rate of V and v is the decay rate of U . The equilibrium 
concentration of b is uo/v^ where uq is the feed rate con- 
stant. The chemical species U and V can diffuse with in- 
dependent diffusion constants and Dy . All the model 
parameters are positive. 

The master equation associated with Eq.(IT]) can be 
mapped to a second-quantized description following a 
procedure developed by Doi p^ . Briefly, we intro- 
duce annihilation and creation operators and a\ for 
V and bi and b\ for U at each lattice site i, with 
the commutation relations [ai,a^] = 5ij and [61,^*^] = 



The vacuum state satisfies ai|0) 



6.10) 



0. We then define the time-dependent state vec- 
tor \^{t)) = Ew,M^({'^>'M,t)n.(«Ir'(6l)"'|0). 
P{{m},{n},t) is the probability to find rm U and Ui 
V particles at site i at time t. The master equation 
can be written as a Schrodinger-like equation — ^^'^[^'^^ = 
H\^(t)), where the lattice hamiltonian or time-evolution 
operator is a hmction of ai^a\,bi, b\ and is given by 
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This has the formal solution |*(t)) = exp(-i?i)|^'(0)). 

The operator Eq. ^ is next mapped onto a continuum 
field theory. This procedure is now standard and we refer 
to [l^ for further details . For the GS system, we obtain 
the path integral exp{~Ht) = / VaVaDbVb e-^^"-''-^-^\ 
over the continuous (and generally complex) stochastic 
fields a{x,t),d{x,t),b{x,t),b(x,t) where the action 5 is 
given by 

S ^ j d"^^ j dt[ddta + DyWdWa + bdtb + DuVWb 
+ fi{d - l)a + v{b - 1)5 - Uo(5 - 1) 
- ^(a-'^a^fo-a^a^M)]. (3) 



We have omitted terms related to the initial state. 
Apart from taking the continuum limit, the derivation 
of this action is exact, and in particular, no assump- 
tions regarding the precise form of the noise are required. 
For the final step, we perform the shift d = 1 + a* 
and 6 = 1 + 6* on the action 5. We represent the 
terms quadratic in a* , 6* by an integration over Gaus- 
sian noise terms, which allows us to then integrate out 
the conjugate fields jl^. To proceed, we note that 
^(XaH(a'^~a*h*)) ^ J V^Vr] , Tj) e'^"' ^+''' where the 
noise functions ^, rj arc distributed according to a double 

Gaussian as P{^, rj) — exp ( — (^, 'q)V^^ ( ^ ) ' ^^^^^ ^ 
the matrix of noise-noise correlation functions 
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Integrating out the conjugate fields a* and b* from the 
functional integral for this shifted action then leads to 
the pair of coupled nonlinear Langevin equations. We 
define the dimensionless fields, u = —b and v = —a, 
dimensionless time, t ~ fit, and spatial coordinates 

, d. The equations describing 
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X'l^ 1 — 1, 

the field dynamics read: 



drV{xi,T) 
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" V^u - w -I- u^u -I- , r) (5) 



drU{xi,T) — V u u—v u-\ \-ri{Xi,T) 



with noise correlations 

{^{xi,T)r]{x„T)) = 
{■q{x„T)r]{x[,T')) = 



(?7(i„r)) =0 

e v^u 5'^{xi ~~ a;^)(5(T — r') 

-| v^u S{xi ~ i-)(5(T - r') 
0. (6) 



and noise strength e = \/A' 



,(d-l)/2 



d — 2 dimensions, e — 



■ III particular, for 

The non-dimensionalization 
of the problem allows us to characterize explicitly the de- 
pendence on the agents diffusivities, reaction rates and 
spatial dimension. For a fixed reaction rate, as we will see 
below, the transition from perfect to the imperfect mix- 
ing regime is associated with an increase in e. Notice that 
77 has zero autocorrelation but non-zero cross-correlation 
with 1^, indicating that this noise is complex. Using the 
Cholesky decomposition, we can express these complex 
noise components as a linear combination of two uncor- 
related (real) white Gaussian noises 9i, 62 thus: 



S.ixi;T) = vy/eu 9i{xi;T) 



(7) 



rj{xi]T) = t) + « 2"^^ ^2(^1; t).(8) 

Thus u and v are also complex fields. Through this pro- 
cedure we can separate the real and imaginary parts of 
the noise and solve numerically the stochastic non-linear 
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reaction diffusion equations ([5]) for the real and imagi- 
nary parts of the fields, [15i] . Now we can obtain numeri- 
cal information from these complex densities. We expect 
the imaginary parts of these fields to be zero on average, 
since the stochastic averages (u) and (v) correspond to 
the physical reactant densities, [16l |. 

In the mean-field approximation, this system possesses 
three homogeneous solutions: one absorbing state R ~ 
{u — ^^,v — 0) and two non-trivial states B± = {u ~ 

— 2 ^ y = ±±] In our simulations wc con- 

sider the d = 2 case with A = 1, D^/Du — 0.5, Uo = v 
and thus the homogeneous solutions will be i? = (m = 
i,« = 0) and B± = (u = i±Vi-VV^ i ^ i i)^ 

state i?+ is globally unstable. The trivial state R is lin- 
early stable and globally attracting for all > and 
IJL> V. The state is stable if > v. In a narrow re- 
gion of parameter space (/i, v) close to — the trivial 
absorbing state R loses stability through a Hopf bifurca- 
tion. In the noise-free case, it is in the vicinity of this 
region where one can find a great variety of spatiotem- 
poral patterns in response to a localized initial spatial 
perturbation [l7[, although no patterns are found if the 
initial condition is homogeneous. Here we will consider 
the evolution of the system described by Eqs. ([5]) with 
initial homogeneous condition and subject to the inter- 
nal noise induced by the reaction, Eqs. ([7]) and ([5]), in 
different mixing regimes. 

In Fig. [T] we display some "snapshots" of the time 
evolution of the real part of the nutrient field u(x,y,t) 
for two different parameter sets and two diffusion rates. 
When displayed in color, blue represents a concentration 

between — and — , green close to — , yellow repre- 
/i /I fj. 

sents an intermediate concentration of roughly ^ and 
red close to i. On a grayscale, lighter grays correspond 
to low concentration and darker ones to high concentra- 
tion. FiglT]-(a) and (b) represent the time evolution of 
the same initial homogeneous condition {u = jj^,v — ^) 
for ij, = 0.095 and v — 0.03. In case (a), a regime of high 
diffusion rates (e = 0.01225), initially there are some 
local spatial fluctuations about the mean value. Then, 
the fluctuations are homogenized due to the high diffu- 
sivity of the reactants and the system evolves towards 
the inactive homogeneous state R of the mean-field ap- 
proach. This is the result we may expect in a perfect 
mixing regime. In case (b), we have increased the noise 
strength (e — 0.0225) by reducing the diffusion rates, i.e. 
to explore deviations from the mean-field results we al- 
low for the imperfect mixing effects: the strong spatial 
incoherent fluctuations give rise to spatial patterns with 
self-replicating and moving globules. In the interior of 
each of these units, in blue, there is a region with sus- 
tained autocatalytic production of v which is causing the 
local depletion of the substrate u. A similar experiment 
is shown in Fig. [T]-(c) and (d), with a new initial condi- 
tion (m = ^,'>^ = ^^)i and parameter set fi = 0.11025, 
!/ = 0.05. The initial condition is again the same in both 



cases, but in the perfect mixing regime where the system 
has a high diffusion rate (and a low noise intensity such 
as e = 0.01) it evolves to the uniform stable active state 
_B_ (sec case (c)), whereas in the imperfect mixing regime 
with low diffusion rates (and higher noise intensities such 
as e = 0.023) the system evolves to a new active pattern 
with globular structures. Thus, in cases (b) and (d) a 
low diffusion rate has induced fluctuations that drive the 
system away from the absorbing state R or the uniform 
blue state _B_ and produces spatial compartmentaliza- 
tion. Due to the imperfect mixing effects the averaged 



FIG. 1: (Color online) Each row of pictures shows -from left 
to right- the time evolution of 5Ru(i, y, r) for a different sim- 
ulation. The initial condition, at r = 0, is shown on the left 
most picture, see text for details, (a) and (b) have the same 
initial condition and reaction parameters, (a) If e ^ 0, which 
corresponds to the perfect mixing regime, the system evolves 
to the uniform absorbing state -R (here shown at r = 50). 
(b) As e increases, in the imperfect mixing regime, the system 
evolves to a new active state with globular replicating struc- 
tures (r — 50). Equivalently, for a new initial condition and 
parameter set, (c) evolves in the perfect mixing regime to the 
uniform active state B- (r = 300) , whereas (d) , in the imper- 
fect mixing regime, evolves to a new active state (r = 300). 

output of the reactor (integrated over x and y) may also 
change with the diffusion rate. In Fig. (upper) we 
show, for fi = 0.0605, = 0.02 and starting with the 
initial homogeneous condition {u — — ^^), how 

for the average of the real part of the field, which cor- 
responds to the density, different mixing regimes lead to 
different reactor outputs. For very high diffusion rates 
(and low noise intensities e), the system evolves to the 
mean-field solution ((3?M)m/, {^v)mf) = notice that 
(div) > (3?u). If the diffusivity of reactants is reduced 
the system evolves to a new equilibrium state with the 
reversed ratio (SRw) < (SRw). Furthermore, after the initial 
transient time, the output of the reactor shows oscillatory 
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behavior which increases in amplitude as we deviate from 
the perfect mixing regime. In Fig. (lower) we show for 
the same simulations the average of the imaginary part of 
the field, confirming that the averaged imaginary part is 
zero up to the statistical error. Thus in any system gov- 
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FIG. 2: (Upper) After a transient time, in the perfect mixing 
limit e = 0.5, the average of the real part of the field tends 
to the mean-field value . However, in the imperfect mixing 
regime, e = 2, these averages deviate significantly from the 
mean-field solution and show oscillatory behavior. (Lower) 
The averaged imaginary parts of the fields are always negligi- 
ble since the averages correspond to physical, real, densities. 



erned by a nonlinear rate law, a knowledge of the bulk or 



average concentrations is not sufficient to predict neither 
the average rate of the reaction, nor the final equilibrium 
spatial distribution or the averaged reactor output. One 
must also know the spatial distribution of the reacting 
material. Furthermore, there is a strong dependence on 
the ratio of the reaction to diffusion rate, and the mean- 
field approach will only be valid in the limit where this 
ratio vanishes. The noise parameter e is a function of the 
spatial dimensionality suggesting that the resulting rates 
of reaction, spatial distribution and global averaged den- 
sities of reactants may change from two to three dimen- 
sions and the simple extrapolation of two-dimensional 
models results may not be adequate. Finally, these re- 
sults show that one can thus tune the mixing efficiently 
to control the composition of the output of the reactor 
and second, that for certain ranges of diffusion rates we 
may discover unexplored dynamical ranges where spatial 
organization takes place. One may also appreciate the 
relevance that the mixing process will have in biologi- 
cal systems, such as epidemics propagation or in viral 
dynamics [l^, where spatial segregation may facilitate 
coexistence. 

This general-purpose, consistent mathematical and nu- 
merical approach allows us to explore rigorously the ef- 
fects of reaction induced fluctuations in imperfect mix- 
ing regimes with complex fields and noises. Dealing with 
the full system we observe that there is multiplicative 
noise acting on both fields. If we had ignored the zero- 
autocorrelated complex field rj in Eq.®, then only the 
V field would have had noise and both field and ^ would 
have been real (fsj . This approximation is not correct in 
the general case and the system outcome may be signif- 
icantly different in non-linear systems such as this one, 
which is quite sensitive to small fluctuations. 
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